;******************************************
; plot regions with climate regime change 
;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;********************************************
 begin

 dirPath = ".../NFood_replication/data/" ; change to local path of your directory

 nmon = 12;
 syr = 2006;
 fyr = 2100;
 csyr = "" + syr
 cfyr = "" + fyr
 nyr = fyr - syr + 1
 dmis = 9.96921e+36

 case = (/"RCP45","RCP85"/);
 ncase = dimsizes(case);
 casename= (/"SSP2","SSP5"/)

 ;Eight active crops in CLM5:
 ; 17/18=Temperate Corn/irrigated, 19/20=Spring Wheat/irrigated, 23/24=Temperate Soybean/irrigated, 41/42=Cotton/Cotton_irrigated,
 ; 61/62=Rice/Rice_irrigated, 67/68=Sugarcane/Sugarcane_irrigated, 75/76=Tropical Corn/irrigated, 77/78=Tropical Soybean/irrigated
 ;Winter wheat (21/22) has not been implemented yet. Only Spring/summer wheat
 ;cft_c3_crop = 1; cft_temperate_corn = 3; cft_spring_wheat = 5 ... cft sequence in crop landunit
 pfts = (/17,18, 19,20, 23,24, 41,42, 61,62, 67,68, 75,76, 77,78 /)      ; crop index among all pfts
 cfts = (/17,18, 19,20, 23,24, 41,42, 61,62, 67,68, 75,76, 77,78 /)-15   ; cft data index
 ncft = dimsizes(cfts);

;read dimension
 fdim = addfile(dirPath+"/land_cover/RCP85.1deg.landarea.nc", "r");
 lat = fdim->lat
 lon = fdim->lon
 nlat = dimsizes(lat);
 nlon = dimsizes(lon);
 dims = (/nlat,nlon/);
;print(lon)
 rad    = 4.0*atan(1.0)/180.0
 clat   = cos(lat*rad); used to estimate global average flux
 clat!0 = "lat"
 clat&lat = lat

 ;directly use the area and land fraction variables for regional sum
 garea = fdim->area  ;sqrkm
 lfrac = fdim->landfrac
 dydx = garea*lfrac*1e6 ;; sqrkm to sqrm
 wgt  = dydx  ;wgt[size (nlat, nlon)]

 ;name coordinates for wgt, for limited area sums
 wgt!0   = "lat"
 wgt!1   = "lon"
 wgt&lat =  lat
 wgt&lon =  lon
 ;usage: input q(time, lat, lon) [size: (ntim, nlat, nlon)]
 ;wgt[lat, lon]: two-dimensional array corresponding to the rightmost dimensions of q
 ;output qSum = wgt_areasum2(q, wgt, opt)  ; => qSum(ntim)


;read data 
 dd2010 = new((/ncase,ncft,nlat,nlon/),"float"); Control period 2006-2015 for RCP85, RCP45
 ;dd2010@_FillValue = default_fillvalue("float")
 dd2050 = new((/ncase,ncft,nlat,nlon/),"float");
 dd2075s = new((/ncase,ncft,nlat,nlon/),"float");
 do kk=0,ncase-1
   ;read crop land cover info from fsurfdata
   ;NOTE: these landuse.timeseries files are very large (24 GB each), and can be downloaded from the NCAR svn site: https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/lnd/clm2/surfdata_map/

   if (kk .lt. 1) then; ;only for RCP45
     fsurf := addfile("/projects/NS2345K/yfan/inputdata/landuse.timeseries_0.9x1.25_SSP2-4.5_78pfts_CMIP6_simyr1850-2100_c190102.nc", "r")
   else ; for RCP85/CCT/MSB/SAI, and RCP45_85LU, RCP85_45CO2
     fsurf := addfile("/projects/NS2345K/yfan/inputdata/landuse.timeseries_0.9x1.25_SSP5-8.5_78pfts_CMIP6_simyr1850-2100_c181209.nc", "r")
   end if
   fcrop:= fsurf->PCT_CROP ; total percent crop landunit of a gridcell; years 1850-2100; [time | 251] x [lat | 192] x [lon | 288]
   fcft := fsurf->PCT_CFT  ; percent crop functional type on the crop landunit (% of landunit); [time | 251] x [cft | 64] x [lat | 192] x [lon | 288]
   fcrop0 := fcrop(156:250,:,:)/100  ;fraction of all crop area 2006-2100
   fcft0  := fcft(156:250,:,:,:)/100 ;fraction of each cft
   ;most crops are within lat[-60:60]

   ;read data for each cft
   do jj=0,ncft-1
   
     ;save annual cultivation area per crop [nyr,ncft,nlat,nlon]
     wgt3d := conform_dims(dimsizes(fcrop0),wgt,(/1,2/)) ; conform wgt to 3d: [year | 95] x[lat | 192] x [lon | 288]
     ;save annual crop planted area per cft; wgt is the land area per grid (static), but fcrop0 and fcft0 are dynamic
     crop_area := wgt3d*fcrop0*fcft0(:,cfts(jj),:,:) ; area per crop per year for all lat/lon
     crop_area := tofloat(crop_area/10000*1.0e-6); sqrm to Mha

     dd0 := crop_area(0:9,:,:)   ;CNTL time period: 2006-2015
     dd1 := crop_area(39:48,:,:) ;future time slice 1: 2045-2054 
     ;dd2 := crop_area(84:93,:,:) ;future time slice 2: 2090-2099
     dd2 := crop_area(69:93,:,:) ;future time slice 2: 2075-2099
     dd2010(kk,jj,:,:) = dim_avg_n_Wrap(dd0,0) ;
     dd2050(kk,jj,:,:) = dim_avg_n_Wrap(dd1,0) ;10-yr avg along the first dim
     dd2075s(kk,jj,:,:) = dim_avg_n_Wrap(dd2,0) ; 
   end do
 end do
 ;printMinMax(dd2075s(0,:,:,:),0)
 ;printMinMax(dd2075s,0) ;[ncase,ncft,nlat,nlon]
 
 ;crop area at the end of 21 century
 ;irrigated vs. rainfed for each crop ; (distinguish irrigated in the global map by slash pattern if >50% of this crop in a grid cell)
 wrk1 = new((/ncase,12,nlat,nlon/),"float"); ;6x2 CPFs irrigated vs. rainfed
 wrk2 = new((/ncase,6,nlat,nlon/),"float");  ;6 CPFs (sum irrigated + rainfed)
 wrk3 = new((/ncase,6,nlat,nlon/),"float");  ;percent irrigated

 ;sum temperate + tropical for corn (0/1) and soybean (4/5)
 wrk1 = dd2075s(:,0:11,:,:) ;excl. tropical corn and soybean, add below
 wrk1(:,0,:,:) = dd2075s(:,0,:,:)+dd2075s(:,12,:,:) ;Rainfed corn (tropical+temperate)
 wrk1(:,1,:,:) = dd2075s(:,1,:,:)+dd2075s(:,13,:,:) ;irrigated corn
 wrk1(:,4,:,:) = dd2075s(:,4,:,:)+dd2075s(:,14,:,:) ;Rainfed soybean
 wrk1(:,5,:,:) = dd2075s(:,5,:,:)+dd2075s(:,15,:,:) ;irrigated soybean

 ;sum irrigated + rainfed for each crop (6 types) ;
 wrk2(:,:,:,:) = wrk1(:,(/0,2,4,6,8,10/),:,:) + wrk1(:,(/1,3,5,7,9,11/),:,:)
 ;percentage of irrigated (fill with pattern)
 wrk2@_FillValue = 0 ;to avoid division by 0 for wrk3, set 0 values of wrk2 to nondata
 wrk3(:,:,:,:) = wrk1(:,(/1,3,5,7,9,11/),:,:)/wrk2
 wrk2@_FillValue = default_fillvalue("float") ;set back for drawing the map
 wrk3@_FillValue = default_fillvalue("float")
  printMinMax(wrk2,0)
  printMinMax(wrk3,0)

 cftname := (/ "Corn", "Spring Wheat", "Soybean", "Cotton", "Rice", "Sugarcane" /)
 ncft2 = dimsizes(cftname);
 cropname = (/"Maize", "Wheat", "Soy", "Cotton", "Rice", "Sugarcane"/)

;##############################
;pick RCP45/RCP85 to plot 6 crops
 casen=1 ;choose a case
 wrk22=wrk2(casen,:,:,:)*1000 ;Mha to kha
 wrk22!0="case"
 wrk22!1="lat"
 wrk22&lat = lat
 wrk22!2="lon"
 wrk22&lon = lon
 wrk22 := lonFlip(wrk22); 0~360 to -180~180
; copy_VarCoords(dd, wrk2);
 wrk33:=wrk3(casen,:,:,:) ;specify the case index
 wrk33!0="case"
 wrk33!1="lat"
 wrk33&lat = lat
 wrk33!2="lon"
 wrk33&lon = lon
 wrk33 := lonFlip(wrk33); 0~360 to -180~180



;plot
 imagePath = dirPath+"/../figures/supplementary"
 imagename = imagePath+"/Fig.S7-8.CropArea_kha_"+casename(casen)+"_irrigation-final"
 wks = gsn_open_wks("pdf",imagename)
 ;gsn_define_colormap(wks, "ncl_default")             ;choose a color palette

 plot  = new(ncft2,graphic) ;crop area; base contour map 
 plot2 = new(ncft2,graphic) ;irrigation; overlay with fill pattern

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False
 res@lbLabelBarOn         = False                 ;turn off individual label bars 

 res@cnLineLabelsOn       = False                 ; turn off line labels
 res@cnFillOn             = True                  ; turn on color fill
 ;res@cnFillPalette        = cmap                 ; cnFillPalette spans the color map by default
 ;res@gsnSpreadColors      = True                 ; this based on an old color scheme that was deprecated
 res@lbLabelAutoStride    = True                  ; optimal labels
 res@cnLinesOn            = False

 res@gsnStringFontHeightF  = 0.020
 res@tmXBLabelFontHeightF  = 0.020
 res@tmYLLabelFontHeightF  = 0.020

; res@cnFillMode = "RasterFill"

 res@mpCenterLonF = 0.
 res@mpMinLatF = -60.
 res@mpMaxLatF =  60.

 res@mpShapeMode  = "FreeAspect"
 res@vpWidthF      = 0.8   
 res@vpHeightF     = 0.4   

 res@tmXBMode = "Explicit"
 res@tmXBValues = (/-180., -120., -60., 0., 60., 120., 180./)
 res@tmXBLabels = (/"180W","120W","60W", "0", "60E","120E","180E"/)
 res@tmYLMode = "Explicit"
 res@tmYLValues = (/-60., -30., 0., 30., 60., 90./)
 res@tmYLLabels = (/"60S","30S","0","30N","60N","90N"/)


; to have a common label bar, all plots should be set to the same interval
; b/c the label bar is drawn by default from the interval of the first plot.
;  res@cnLevelSelectionMode =  "ManualLevels"   
;  res@cnMinLevelValF       = 0  ;Mha
;  res@cnMaxLevelValF       = 1 
;  res@cnLevelSpacingF      = 0.1
; difficult to assign white color to small values>0 without using ExplicitLevels. 
; change a default color map
;  ;---Reading the default colormap into a two dimensional 254 x 4 (RGBA) array
;  ;cmap = read_colormap_file("rainbow")
;  cmap = read_colormap_file("ncl_default")
;  ;cmap = cmap(::-1,:) ; reverse color map
;  ;select the middle cmap elements that correspond to 0 in the middle of value range (-0.5~0.5)
;  ; cmap(120:136,:) = 1 ; use white for zeros at the middle; 0 for black, 1 for white
;  cmap(0:25,:) = 1 ; use white for zeros at the begining; 0 for black, 1 for white
;  res@cnFillPalette         = cmap     ; use the modified color map; cnFillPalette spans the color map by default
  
 ;##########colors for crop area
 res@cnLevelSelectionMode = "ExplicitLevels"
 res@cnLevels             = array_append_record((/1/), ispan(50,500,50), 0) ; 11 level to divide data into 12 ranges
 res@cnFillPalette        = read_colormap_file("precip_11lev") ;12 colors, first white
 res@cnMissingValFillColor = "white"  ; for missing values 0 (default gray)
 ;res@cnFillColors = (/ 0, 54, 2, 3, 4, 5, 28, 20, 14, 0/)
 ;res@cnSpanFillPalette     = True    ; spans cnFillPalette to get color index values for cnFillColors in Explicit mode

 ;#########fill patterns for irrigated area
  res2                      = True              ; set up a second resource list
  res2@gsnDraw              = False             ; do not draw the plot
  res2@gsnFrame             = False             ; do not advance the frame

  res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res2@cnMinLevelValF      = 0.0         ; set min contour level
  res2@cnMaxLevelValF      = 1.0         ; set max contour level
  res2@cnLevelSpacingF     = 0.1         ; set contour spacing

  res2@cnInfoLabelOn       = False       ; turn off info label
  res2@cnLinesOn           = False       ; do not draw contour lines
  res2@cnLineLabelsOn      = False       ; do not draw contour labels
  
  res2@cnFillScaleF        = 0.5         ; add extra density
  res2@cnFillOn            = False
  res2@cnFillColor         = "red"       ; a get around way to specify a pattern fill color into gsn_contour_shade.

;############ Overlaying a stipple/point pattern to show area of interedst
;adapted from https://www.ncl.ucar.edu/Applications/Scripts/conOncon_4.ncl
; or more fency plot: http://www.ncl.ucar.edu/Applications/Images/coneff_14_lg.png
;############

 do kk=0,ncft2-1
   res@gsnCenterString     = cropname(kk)
   plot(kk)                = gsn_csm_contour_map(wks,wrk22(kk,:,:),res)      
   plot2(kk)               = gsn_csm_contour(wks,gsn_add_cyclic_point(wrk33(kk,:,:)),res2)
   opt     = True
   opt@gsnShadeFillType    = "pattern"
   ;opt@gsnShadeMid         = 17  ;points ; pattern fill areas >= lowval contour level and <= highval
   opt@gsnShadeHigh        = 6  ;pattern #6 cross hatch 45 angle; see NhlTFillIndex 
   opt@gsnShadeFillScaleF  = 0.6 ;a value <1 will make the fill pattern  denser
   plot2(kk)   = gsn_contour_shade(plot2(kk), 0.1, 0.5, opt) ;shade contour levels >=highval 0.5 for irrigation
   
   overlay(plot(kk),plot2(kk))  ;overlay the two plots   

   ;plot(kk) = gsn_csm_contour_map_overlay(wks,wrk22(kk,:,:),\
   ;           wrk33(kk,:,:),res,res2) ;draw two maps together
   ;draw(plot(kk))           
 end do


;Draw Labelbar for whole panel
;*********************

;panel
 panres1                  = True
 panres1@gsnFrame         = False
 panres1@gsnPanelBottom   = 0.05          ;add some pace in the bottom
 panres1@gsnPanelYWhiteSpacePercent = 5   ;add 5% white space between plots
 panres1@gsnPanelMainString = casename(casen)+" crop area by the end of 21st century (2075-2099)"     ; set panel title
 panres1@gsnPanelLabelBar = True ;add a common label to whole panel
 panres1@lbLabelFontHeightF = 0.011 ;set label tick font size
 gsn_panel(wks,plot(0:5),(/3,2/),panres1)        ; now gsnDraw is on

; Draw a text string at the bottom below label bar
 txres               = True
 txres@txFontHeightF = 0.015
 gsn_text_ndc(wks,"Crop Area (x10~S~3~N~ ha)",0.5,0.05,txres) 
  ;function codes: ~F8~D~ is for big delta; ~F~ means switch back to default font
  ;gsn_text_ndc(wks,"~F8~D~",0.5,0.1,txres) ;~F8~D~ is the special code for delta

 frame(wks)


end
